COMMENT ⊗   VALID 00005 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	C This is taken from the "Handbook of Mathematical Functions"
C00005 00003		SUBROUTINE SIDEBANDS(CAR,XMOD,CINDEX,K)
C00006 00004		SUBROUTINE REFLECT(MAX,K)
C00008 00005
C00010 ENDMK
C⊗;
C This is taken from the "Handbook of Mathematical Functions"
C by Abramowitz and Stegun, Dover press S1272, 1965.
C Chapter 9 is on Bessel Functions. The polynomial
C approximations are found in section 9.4, pages 369
C and 370. The recursion formula for generating higher
C orders from J0 and J1 is found at the bottom of
C table 9.1, page 390

	SUBROUTINE BESCOEF(MI,X,I,K)
	REAL MI,J0,J1
	COMMON FREQ(3,0/50,50),FUNC(50)
	IF(I.GT.1)GO TO 30
	IF(MI.GE.3.0)GO TO 10
	xs=(MI/3)**2
	j0=1.0+xs*(-2.2499997+xs*(1.2656208+xs*(-0.3163866+
	1xs*(0.0444479+xs*(-0.0039444+xs*0.00021)))))
	j1=MI*(.5+xs*(-0.56249985+xs*(0.21093573+
	1xs*(-0.03954289+xs*(0.00443319+xs*(-0.00031761+
	2xs*0.00001109))))))
	GO TO 20
10	xs=3.0/MI
	j0=(1.0/sqrt(MI))*(0.79788456+xs*(-0.00000077+
	1xs*(-0.0055274+xs*(-0.00009512+xs*(0.00137237+
	2xs*(-0.00072805+xs*0.00014476))))))*
	3cos(MI-0.78539816+xs*(-0.04166397+
	4xs*(-0.00003954+xs*(0.00262573+
	5xs*(-0.00054125+xs*(-0.00029333+
	6xs*0.00013558))))))
	j1=(1.0/sqrt(MI))*(0.79788456+xs*(0.00000156+
	1xs*(0.01659667+
	1xs*(0.00017105+xs*(-0.00249511+xs*(0.00113653-
	2xs*0.00020033))))))*
 	3cos(MI-2.35619449+xs*(0.12499612+xs*(0.0000565+
	4xs*(-0.00637879+xs*(0.00074348+xs*(0.00079824-
	5xs*0.00029166))))))
20	FREQ(2,0,K)=J0
	FREQ(2,1,K)=J1
	FREQ(2,2,K)=J1
	RETURN
30	Y=I
	IF(MI.LE.0.0000001)GO TO 50
	IJ=I*2		
	X=((2.0*(Y-1.0))/MI)*FREQ(2,IJ-2,K)-FREQ(2,IJ-4,K)
	RETURN
50	X=0
	RETURN
	END
	SUBROUTINE SIDEBANDS(CAR,XMOD,CINDEX,K)
	COMMON FREQ(3,0/50,50),FUNC(50)
	IT=-1.0
	MAX=CINDEX+5.0
	IF(FREQ(1,50,1).LT.MAX*2.)FREQ(1,50,1)=MAX*2.
	DO 100 J=0,MAX*2,2
	I=J/2
	NORDER=I
	IF(I.NE.1)CALL BESCOEF(CINDEX,COEF,NORDER,K)
	IF(I.GT.0)GO TO 20
	FREQ(1,I,K)=CAR
	FREQ(3,I,K)=0.0
	GO TO 30
20	XI=I
	YMOD=XI*XMOD
	FREQ(1,J-1,K)=CAR+YMOD
	FREQ(1,J,K)=CAR-YMOD
	IF(I.EQ.1)GO TO 10
	FREQ(2,J-1,K)=COEF
	FREQ(2,J,K)=COEF
10	FREQ(3,J-1,K)=I
	IF(IT)I=-I
	FREQ(3,J,K)=I
	IT=-IT
30	CONTINUE
100	CONTINUE
102	FORMAT(3F)
	NAX=MAX
	CALL REFLECT(NAX,K)
	RETURN
	END

	SUBROUTINE REFLECT(MAX,K)
	COMMON FREQ(3,0/50,50),FUNC(50)
	DO 100 N=0,MAX*2
100	IF(FREQ(3,N,K).LT.0.0)FREQ(2,N,K)=-FREQ(2,N,K)
	DO 201 J=0,(MAX*2)-1
	DO 201 I=J+1,MAX*2
	IF(FREQ(1,J,K).OR.FREQ(1,I,K).EQ.99999.)GO TO 201
	IF(ABS(FREQ(1,J,K)).NE.ABS(FREQ(1,I,K)))GO TO 201
	IF(FREQ(1,J,K).GT.FREQ(1,I,K))GO TO 20
	FREQ(2,I,K)=FREQ(2,I,K)-FREQ(2,J,K)
	FREQ(1,J,K)=99999.
	GO TO 201
20	FREQ(2,J,K)=FREQ(2,J,K)-FREQ(2,I,K)
	FREQ(1,I,K)=99999.
201	CONTINUE
	RETURN
	END


C     	MAIN PROGRAM
	COMMON FREQ(3,0/50,50),FUNC(50)
10	ACCEPT 100,C,XM,ZI1,ZI2
100	FORMAT(4F)
	CALL GEN
	DO 300 IS=1,25
	CI=ZI1+(ZI2-ZI1)*FUNC(IS)
300	CALL SIDEBANDS(C,XM,CI,IS)
	CALL DISP
	MAX=FREQ(1,50,1)
C	TYPE 101,(((FREQ(K,L,M),K=1,3),L=0,MAX),M=1,25)
101	FORMAT(1H (3F/))
	GO TO 10
	END

	SUBROUTINE GEN
	COMMON FREQ(3,0/50,50),FUNC(50)
	X=1./25.
	DO 100 I=1,25
	Y=I-1
100	FUNC(I)=X*Y
	RETURN
	END

	SUBROUTINE DISP
	DIMENSION I(1),IJ(1000)
	COMMON FREQ(3,0/50,50),FUNC(50)
	CALL DPYSET(1,IJ,1000)
	CALL ALINE(-400,0,400,0)
	CALL ALINE(-400,400,-400,0)
	MAX=FREQ(1,50,1)
	KL=1
	DO 200 J=0,MAX
	IF(FREQ(1,J,KL).EQ.99999.)GO TO 200
	IX=ABS(FREQ(1,J,KL))-400.
	ZZ=IX
	IY=(ZZ+400.)*(-1.)+400.*FREQ(2,J,KL)
	BASE=(ZZ+400.)*(-1.)
	KL=KL+1
	CALL AIVECT(IX,IY)
	DO 200 K=2,25
	IF(FREQ(1,J,K).EQ.99999.)GO TO 200
	IF(FREQ(1,J,K).EQ.0.0)GO TO 200
	IX=IX+20
	IY=FREQ(2,J,K)*400.+BASE
	IBASE=BASE
	IF(IY.LE.IBASE)IY=(IABS(IY)-IABS(IBASE))+IBASE
	CALL AVECT(IX,IY)
200	CONTINUE
	CALL DPYOUT(1)
	RETURN
	END